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AN OPTIMAL COMPACT STORAGE SCHEME FOR 
NONLINEAR REACTOR PROBLEMS BY FEM 



by 

D. Salinas, D. Nguyen and R. Franke 

Abstract - This work shows that optimal compact storage of coefficient 
matrices affords a significant reduction in core storage requirements over 
banded storage schemes. The resulting savings enables in-core finite 
element solutions of large systems not otherwise possible. It is shown 
that Gears method for the stiff system of a nonlinear reactor dynamics 
problem is not as efficient as Crank-Nicolson integration because of 
substantially greater core requirement, despite its superior tracking 
ability. A remedy in the form of a modified implicit version of Gears 
method with a significant reduction in core requirements is shown to pro- 
vide the same excellent accuracy as Gears method. Comparisons between the 
modified Gear method and the Crank-Nicolson method show the relative 
advantages and disadvantages of each. Finally, it is shown that although 
the nonlinearity encountered in this problem can be treated directly, a 
linear approximation of the nonlinear term affords a substantial reduc- 
tion in core requirement with a relatively small cost in accuracy. 

Statement of the Problem 

The finite element method (FEM) has been successfully applied to the 
problem of a space-time nonlinear reactor with temperature dependent feed- 
back (1,2). However, the effective numerical treatment of the resulting 
system of stiff O.D.E. called for further investigation. This work com- 
pares alternative computational schemes in a FEM solution of a nonlinear 
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nuclear reactor problem. 

A test case is chosen which consists of a fully-reflected fast 
reactor, subject to super prompt critical conditions. By neglecting the 
delayed neutrons and using the prompt feedback model, the following non- 
linear reactor dynamics equation is obtained in (r,z) space in the core, 

1 || (r,z,t) _ Da2 ^ + XE ^ . wtp 2 = 0 (i) 

A similar equation without the nonlinear feedback term results in the 
reflector. A detailed development of these equations is given in ( 1_) . 

In accordance with Galerkin's method, Eq. (1) is transformed into 
the implicit system of N ordinary differential equations, 

£ - B^tj) + £ £ C iJk *j* k - 0 (2) 

j=l j=l k=l 

i , = 1 , . . . , N 

together with appropriate initial conditions, <|^(0), i = 1, ..., N. The 
computational efficiency of any finite element program depends on the al- 
ternative schemes available for selection by the investigator. The pre- 
sent work considers the relative merits of various integration methods, 
algebraic equation solvers, array storage methods, as well as the effects 
of linearization. 

Linearization 

A well known procedure for treating nonlinear terms in FEM is by a 
"linearization" of the nonlinear terms. Linearization here refers to the 
process whereby a nonlinear term F [^(t)] is approximated by F (t-At) ] , 
where <j;. (t-At)is the known value of at the previous time. For the par- 
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N N 

icular quadratic nonlinear term z z C. -•■,ij>. i (t)i|»| / (t) of Eq. (2), a 

j=l k=l 1JK J k 

better linear approximation is 

N 

z C* . \p -(t) , where 
j=l 1J J 

C i j = C i jk^k^ ^ 

It should be noted that although the nonlinear quadratic term ssC -j jk^j^k 
can be treated directly, and indeed it was in reference (1_) , the linear 
approximation, Eq. (3), results in a substantial reduction in both array 
storage and CPU time for a relatively small expense in accuracy, less than 
5 percent for the cases run. The savings in storage and CPU time depends 
on the integration method and the algebraic equation solver employed, as 
well as the method of array storage. 

By far the most common (but not optimal) method of array storage is 
the band storage scheme, whereby only the bandwidth (or half bandwidth in 
the case of symmetric matrices) terms in a matrix are stored. Since ma- 
trices obtained in FEM analyses are always banded, banded storage is always 
advantageous. Table 1 lists the storage requirements and the storage 
savings accomplished by linearization of symmetric and nonsymmetric C^.^. 
Coefficient a is the number of bytes per word for the particular computer 
being utilized. For an IBM360 computer, a=4 for single precision, and 8 
for double precision. N is the number of unknowns and n is the total 
bandwidth. 
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TABLE 1. 



Matrix 

Condition 


C ijk 




c*. 

ij 


1 

Reduction 

1 


non sym.,non banded 


aN 3 


i 


aN 2 


aN 2 (N-l) 


non sym. , banded 


aNn 2 


i 


aNn 


aNn(n-l ) 


sym.,non banded 


2f<N+l)(I + 


2N+1 \ 
6 ‘ 


aN(N+l ) 
2 


a(N+l )N /2N+1 
2 ' 6 V 


sym. , banded 


i^jkn+l ) + 


2n+l \ 
6 } 


2f(n*l) ! 


a N/__ Ll \ /2n+l 1\ ’ 

-T n+1 H 6 ‘ 2 } 



It should be pointed out that some additional reduction can be obtained 
for symmetric matrices by discarding the n - ^ — zeroes in the lower right 
triangle of the banded N x n Ct. matrix. In most cases n is small relative 
to N and the - n - Q>~ -- ) - additional reduction is not worth the effort. 

Table 2 gives an idea of the savings achieved using banded storage 
for some typical regular rectangular grids as shown in Figure 1. The 
results are associated with the symmetric C — ^ and matrices of the non- 
linear reactor dynamics problem considered in this work, in which linear 
triangular elements were utilized with single precision on an IBM360 com- 
puter. (R-l) is the number of rows and (S-l) the number of columns in the 
rectangular grid, with S > R. 
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TABLE 2. 



Grid" 

• s 


c. .. 

ijk 


c*. 

ij 


Reduction 


1 R=5 S=8 

! N=40 n=3 
1 


5,600 


1,120 


4,480 


! R=10 S=20 
j N=200 n=23 


80,000 

I 


9,600 


70,400 


R=20 S=40 
N=800 n=43 


1 

1 

1 ,056,000 


70,400 


985,600 


R=40 S=60 

;N=2400 n=83 


11 ,424,000 


403,200 


11 ,020,800 



Optimum Compact Storage of Coefficient Matrices 

J.L. 

Since the interpolation function for the i L unknown is non-zero only 
over those elements which contain i (i.e., the basic functions in FEM 
have local support only), the matrices resulting from FEM are both banded 
and sparse. The sparce nature arises from the fact that the unknowns 

i.L 

surrounding the i tn unknown are not consecutively numbered. Common prac- 
tice takes advantage of the banded nature of these matrices via a banded 
storage scheme, which stores an NxN matrix as an Nxn matrix where N is 
the number of unknowns and n is the total bandwidth of the matrix. 

Using banded storage, special care in the numbering of unknowns can affect 
a substantial decrease in n. As a simple example of this, consider the 
regular rectangular grid of Figure 1, having (R-l) rows of elements, and 
(S-l) columns of elements. For linear triangular elements, sequential 
vertical numbering results in n = 2R + 3, whereas sequential horizontal 
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numbering gives n = 2S + 3. The saving achieved by vertical numbering, 
assuming R < S, is 2N(S-R) words, or 2aN(S-R) bytes, where a is the num- 
ber of bytes per word. For symmetric matrices the previous numbers are 
hal ved. 

These results may be generalized for triangular elements of order t. 
For the rectangular grid of Figure 1, the bandwidth for triangular ele- 

j_ i_ 

ments with t in order polynomial interpolation is 

n = 6Rt - 4R - 2t + 5 (4) 

for non symmetric matrices, and the half bandwidth is 

r s = 3Rt - 2R - t + 3 (5) 

for symmetric matrices. Equations (4) and (5) are derived on the basis of 
condensing out the interior element nodes. Table 3 presents the evalua- 
tions of these formulas for linear, quadratic and cubic elements for the 
rectangular grid of Figure 1. 

TABLE 3. 



t 


n 


n 






s 


1 


2R + 3 


R + 2 


2 

i 


8R + 1 

j 


4R + 1 

i 


3 


14R - 1 

; 


7R 
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For banded systems of linear algebraic equations, direct methods of 
the Gaussian elimination type are most often utilized, being preferred over 
iterative methods such as successive over-relaxation (SOR). Depending on 
the particular application, there are points in favor of each method. In 
the case of linear equilibrium (elliptic) boundary value problems Gaussian 
elimination has the following advantages: 

1. Direct methods allow the use of iterative refinement. 

2. If more than one right hand side must be processed, the initial 
cost of decomposition must be expended only for the first solu- 
tion. 

For transient (parabolic or hyperbolic) initial-boundary value problems, 
or multi-step nonlinear equilibrium problems, SOR has the following 
advantages: 

1. The large number of iterations for a solution at the first step 
is expended only once, succeeding solutions require far fewer 
iterations since a good estimate of the solution is then avail- 
able. 

2. No working space is required by SOR and therefore matrices may 
be compacted to maximum density. 

It is this last point which suggests the replacement of the banded storage 
schemes by optimum compact storage (OCS) which will be described shortly. 
The argument in defense of direct methods that storage is becoming in- 
creasingly abundant is a spurious one which impedes investigation of 
larger systems, especially for three dimensional problems. It will be 
shown here that optimum compact storage of coefficient matrices leads to 
a substantial saving of core for large systems. A description of OCS 
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follows . 



Accounts of OCS are given by Tewarson (_ 3 ) and Gustavson ( 4 J . The 
underlying idea behind OCS is simply to store only the non zero terms of 
a coefficient matrix. The reduction in storage requirement achieved 
through OCS as compared to banded storage is illustrated by the regular 
rect angular grid shown in Figure 1. It is striking that regardless of 
R and S, the band of n terms surrounding any nodal point contains, at most, 
7 non zero coefficients in the case of linear triangular elements. For 
the cases of quadratic and cubic elements the maximum number of non zero 
coefficients for any nodal point is 19 and 30 respectively. Hence the 
maximum storage requirements for the regular grid of Figure l for linear, 
quadratic and cubic triangular elements are a(7N), a(19N) and a(30N), 
respectively. The actual reduction is not simply the difference between 
banded core requirement and OCS core requirement for two reasons. 

1. OCS storage can be achieved with less than cc(yN) where y is 
the maximum number of 'non zero coefficients and a is the number 
of bytes per word. 

2. Additional vector arrays are required for the implementation of 
OCS. 



Although the following discussion is given in terms of the linear 

NxN C*. • array, the nonlinear NxNxN C • .. array is treated in the same way. 
1 J l J K 

Indeed the direct treatment of the nonlinearity by OCS provided the most 
impressive reduction in core imaginable. An additional word on this will 
be given at the end of this section. It is well to keep in mind, however, 
that although direct treatment of the nonlinearity is esthetically pleas- 
ing, it is a questionable luxury in view of the excellent results obtained 
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by linearization. 



The implementation of OCS requires two integer array vectors, say 
ISTART and NAME, and a vector of the non zero coefficients of C, which 
we shall call CC. The i^ integer entry in the (N+l)xl ISTART vector 

is the number q . , where q. = ^ p.+l , and p. is the number of terms 

1 1 i=l 1 1 

th J 

in the i equation, i.e., p^ is the number of nodes surrounding the 

i^* 1 node. ISTART, then, is a pointer vector whose i^ term locates the 

initial position in the CC vector of the contributing coefficients to 

the i^* 1 equation. The Mxl NAME vector, where M = ^ p . , is composed of 

i=l 1 

N successive vector blocks of variable length p. , i=l,...,N. The p^ 

i. L, 

integer numbers in the i tn block of NAME identify the p^ contributors to 

J. L. 

the i tn equation. The Mxl CC vector contains the real non zero coeffi- 
cients of the NxN matrix C, arranged in the same contiguous block arrange- 
ment as the NAME vector. The term in the i^ block, CC (ISTART 
(I) + J-l) is coefficient Cj k where K = NAME (ISTART (I) + J-l). To fix 
ideas the vectors associated with the simple grid shown in Figure 1 are 



ISTART = < 1 


5 10 13 18 


25 


30 33 


38 


41 


> 


NAME = < 1 


2 4 5| 2 3 1 6 


5 1 


1 


9 


8 7 


> 


O 

O 

II 

A 

O 


C 12 C 14 C 15 ! C 22 C 23 


C 21 


C 26 C 25 1 * 


• • 


1 C 99 


C 98 C 97 > 



In this case N=9 and M=40. 

With the interior nodes of an element condensed out, the length 
(number of entries) of the NAME and CC vectors for the regular rectangular 
grid of Figure 1 with triangular elements of t Ln order polynomial interpo- 
lation, is: 



9 



(5) 



N 

M = E p = RS ( 1 5 t 2 - 6t - 2) + (R + S) (-14t 2 + 8t + 2) 
i=l 

+ (13t 2 - lOt - 1) 

The number of unknowns is 

N = RS (3t - 2) + (R + S) (2 - 2t) + (t - 1) (6) 

Table 4 gives the evaluation of Equations (5) and (6) for linear, quadratic 
and cubic triangular elements for the grid of Figure 1. 

TABLE 4. 



t 


N 


M 


» ! 
! i i 

S 


| RS 

i 


7RS - 4 (R + S) + 2 


2 


4RS - 2 (R + S) + 1 


46RS - 38 (R + S) + 31 


3 


7RS - 4 (R + S) + 2 


115RS - 120 (R + S) + 86 



Hence the core requirements for banded and compact storage are: 

N n* = n for non sym. matrices 

B 

n* = n/2 for sym. matrices 

(7) 

N * oM + e (M + N + 1) 

respectively, where a and b are the number of bytes per word for real and 

integer numbers. For an IBM360, a is 4 for single precision, and 8 for 

1 5 

double precision; B is 2 for integers less than 2 and 4 for integers 
1 5 

greater than 2 . Comparative results between banded and OCS storage, 
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using the formulas from Table 2 and 3, for single precision with an IBM360 
computer for symmetric operators on some typical RxS rectangular grids are 
given in Table 5. 

TABLE 5. 



r — 

I 

R = 10 S = 20 

l 


I 

R = 20 S = 40 


» 

i 

R = 40 S = 60 


t | 1 

i 


2 


3 


1 


I 

2 ! 3 

: 

. 


i 


2 3 j 


! 

| N i 200 

i ! 


741 


1 ,282 


800 


3,081 j 

1 


5,362 

i 


2,400 


l •; 

9,401 j 16,402 

1 


M ! 1,282 

{ 


8091 


19,486 

i 




5,362 


i 

34,551 j 84,886 

i 


16,402 


106 ,686 i 


264,086 


' 

n i 12 

i s 


41 

; 


7oj 22 

i 


*81 140 


42 


161 


280 ' 

{ 

J 


; n b 


9,600 121 ,524 


358,960; 70,400 


998,244 '3,002,720 

t 


403,200 


6,054,244 


r 

18.4x10° 


* 

N c 


8,094 50,030 


119,482 


33,774 


' p " 

213,470 520,042 


103,214 


658,920 


1 ,617,322 


«B- N c 


1,506 


71 ,492 


239,478 


1 

36,626 


784,774*2,482,678 


299,986 


5,395,324 


1 6 . 8x1 0 6 


°/o Dif 

< 

1 


-15.7 


-58.5 


-66.7 


-52.0 


-78.6 


-82.7 


-74.4 


-89.1 


-91.2 



The notation in Table 5 is as follows: 

t - order of polynomial interpolation 
N - number of unknowns 

M - number of terms in the NAME and CC vectors 
n $ - bandwidth for the symmetric matrix 

Ng - core required for banded storage, single precision IBM360 



11 



N c - core required for OCS, single precision IBM360 

N b -N c - reduction in core storage by OCS, single precision IBM360 

% Diff - percent difference between banded storage and OCS 

In Table 5, there is no intention to suggest comparisons for a par- 
ticular grid with respect to the order of polynomial interpolation. Com- 
parisons with respect to N, the number of unknowns, are more meaningful. 
ISTART and NAME vectors are the same for the A and B matrices and hence 
the comparisons given in Table 5 are conservative for systems with more 
than one coefficient matrix. 

The algorithm to assemble the element matrices into the CC vector 
is straight forward and can be accomplished in a dozen or so steps. In 
our work it was found useful to construct and have the computer punch 
out the ISTART and NAME vectors from the element-system connectivity of 
the discretized model under consideration. Again this pre-processing 
algorithm was quite simple to program. 

Solution Methods 

It was not the intention here to compare the many techniques avail- 
able for the solution of Eq. (2), but rather to establish the computational 
features of the multi-step, predictor-corrector Gear method for stiff 
systems, and the single-step, implicit Crank-Nicol son method. The work 
led to the development of a modified implicit form of Gears method which 
incorporated some of the attractive features of each system. 

Having been specifically developed for stiff systems. Gears method 
handled Eq. (2) quite efficiently with respect to CPU time. Indeed, Gears 
method achieved a specified accuracy criterion with far fewer integration 
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steps than the Crank-Nicolson method, which experienced particular diffi- 
culty during the early transient period of time. However, Crank-Nicolson 
performed equitably with Gears method during the later transient period 
as shown in Figure 2. 

A severe drawback to multi-step methods such as Gears method, in 
contrast to single-step methods such as Crank-Nicol son, is the need to 
retain previous time solutions. The predictor-corrector Gear method, as 

•hh 

typified in the DVOGER subroutine (5.) is capable of 6 tn order interpolation 
in time and requires N(N + 17) words of work storage, in addition to 8N 
words for storage of past history, unquestionably a considerable expense. 

Another disadvantage of Gears method is the need to transform the 
implicit system of Eq. (2) into the explicit system 

\p = A"^ (B - C^i)tp = F(i|>) (8) 

in the case of direct treatment of the nonlinear term, or 

= A -1 (B - C*)u> (9) 

when the nonlinear term is given the linear approximation of Eq. (3). Equa- 
tions (8) and (9) are formal statements since actual pre-multiplication 
by A~^ is not necessary. These equations may be obtained by a Cholesky 
decomposition into 

LL T * = Bif; - Cn> 2 (10) 

or, 

LL T tji = (B - C)H> (11) 

which, in turn, may be solved by forward and backward substitutions for 
iji. It should be noted that Cholesky decomposition eliminates the possible 
use of OCS, and hence a severe penalty in coefficient array storage is 
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incurred for the A matrix. This, in addition to the need to store the 

3F i 

NxN Jacobian matrix -r— gives a great advantage to the Crank-Nicol son 
method, despite the superior ability of Gears method for stiff systems. 
The modified Gears method, to be described shortly, offers a substantial 
dimini shment in the storage requirement of Gears method, while maintain- 
ing the advantage of superior integration for stiff systems. The net 
effect is a method which compares favorably with the Crank-Nicol son 
method. 

The Modified Gear Method 

Our initial work (J_) presented a dilemna which required additional 
investigation. On the one hand the superior tracking ability of Gears 
method for the stiff system of equations (2) was admired, while on the 
other hand there was dissatisfaction with the large core storage require- 
ment. One way to minimize this difficulty was a modification of Gears 
method. 

1. To treat the implicit system of Eq. (2) and thus achieve the 
core saving advantages of OCS, and 

2. Eliminate the need for storing the NxN Jacobian matrix. 

In order to obtain an understanding of the modified Gear method 
the following presentation is given. 

The core storage requirements of the multi-step Gear method arise 
primarily from the necessity to compute and store the NxN Jacobian matrix 
9Fi/8iK . While A, B and C (or C*) are sparse banded matrices, the 

J 

Jacobian A"^ ( B-2Cip ) , or A”^(B-C*), is full. Storage of this matrix 
2 

accounts for N of the work storage words required by DVOGER. At the 
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suggestion of Prof. C. W. Gear (6), an existing program (7) for implicitly 
defined stiff differential equations was modified to take maximum advantage 
of the sparse nature of matrices A, B, and C (or C*) . This includes the 
use of OCS for the equivalent of the Jacobian matrix, which uses the same 
ISTART and NAME vectors as the A and B matrices, and the use of iteration 
(SOR) to solve for the difference between the predicted value and the 
corrected value at each time increment. 

A brief background concerning Gears method for implicit stiff equations 
is necessary in order to understand why SOR results in a substantial increase 
in computational efficiency for problems of this type. In the general case, 
one seeks to obtain a numerical solution of the system of equations 

F (t, y, y) = 0 (12) 

with initial conditions y(t Q ) specified. At each time step here denoted 
by superscript k, the corrector equations 

F (t, y (k) , - V (k) /3 oh + £ k ) = 0 (13) 

♦ (kl 

which result from approximating y by a linear combination of y v ' and 

f kl 

previous solution values, must be solved for y v ' . In the above, h is the 

stepsize, a Q and e Q depend on the order of the method (interpolation on 

time), and depends on the order and previous solution values. Gears 

( kl 

method uses Newtons method to solve for y v , hence the Jacobian, 

ot 

J = -j ^7 - (-r\) ■— , (or a reasonable approximation to it) must be computed. 

°y p n n 3* 

0 y (k) 

Because the predictor provides a good estimate of y v ' and Newton's method 
solves for the difference in iterates, <5y, the linear equations to be 
solved are of the form 

J ( <5y ) = -F. (14) 
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When the solution <5y is added to the predicted value it is apparent that 
6y need be accurate to only a few significant digits. Also, because up 
to three Newton iterations are performed, inaccuracies in the first 
iteration, where the largest corrections, Sy, are expected, are auto- 
matically accounted for in succeeding iterations. It should be noted 
that the iteration used by Gear is actually a modi fed Newton method since 
J is not evaluated at each iteration, nor even at each time increment. 

The error tolerance in solving for Sy by SOR was satisfied when 
either of two criteria were met 

1. The difference in components of successive iterates were as small 

compared to the solution Sy as that allowed by Gears method for 

( kl 

the relative error in y v , or 

2. The difference in components of successive iterates compared 

(k) 

to y v ' were less than 1 percent of the error allowed by Gears 
f kl 

method in y v . 

(kl 

Considering the comparitive magnitudes of the vectors y v ' and 8y, the 
result was that very few SOR iterations are required to solve equation 
(14). In the problem under investigation here, 3-4 iterations were 
required for the first Newton iterate, and 1-2 for second and (if required) 
third Newton iterates. 

A comparison of the number of arithmetic operations required for 
solution of the corrector equation (13) for three different schemes 
follows. An operation is defined to be a multiplication or division 
and an addition or subtraction. 

1. Use of SOR in the modified Gear method results in M + 6N opera- 
tions per iteration of SOR, including those required for 
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convergence tests. An overhead of N divisions for the initial 
estimate of the solution (taken to be 6y^/J^) is required for 
each Newton iteration. Storage for J requires only M locations 
for the coefficients since the ISTART and NAME vectors are the 
same as for the A and B matrices. 

2. Use of Cholesky decomposition to factor J and solution of equa- 
tions (14) by forward and backward substitutions requires 

2Nn - n(n-l) operations for each Newton iteration, plus an 

2 

overhead of about Nn operations for Cholesky decomposition each 
time J is computed. While J is only recomputed every few time- 
steps, for large problems the decomposition of J will probably 
take as many operations as all the succeeding solutions of 
equation (14). The matrix J is assumed to be stored in symmetric 
band storage form, requiring Nn storage locations. 

3. In a routine such as DVOGER where the explicit differential 
equations must be used, the banded nature of the Jacobian is 
destroyed and the resulting J matrix is full and nonsymmetric. 

p 

Gauss elimination for the solution of (14) requires N operations 
for forward and backward substitutions, with the LU decomposition 

3 

requiring about N /3 operations each time J is computed. Since 

2 

J is full, N storage locations are required for it. 

The total amount of work involved in solving the corrector equation 
(13) is difficult to compare in the first two instances because of a 
varying number of iterations required in solving (14) by SOR, as well as 
the variable nature of M compared to Nn. It is clear that DVOGER and 
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similar codes which must have the explicit equation, i.e., (8) or (9), 
to solve, cannot compare favorably with the modified Gear method which 
attacks the implicit equation (2), since both storage and run time re- 
quirements are greater. For the problems under consideration here, typ- 
ical test runs indicated the solution of (14) by SOR is much faster than 
elimination, since the total run time for solution of the differential 
equation was less by a factor of about 2.5. With regard to storage, the 
advantage of SOR over Cholesky decomposition is greater than shown in 
Table 5 since the ISTART and NAME vectors are the same as for the A and 
B matrices. The details are presented in ( 8 ). 
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Figure 1. Rectangular Domain with Triangular 
Elements (S > R) . 




Figure 2. Plot of Center Point Neutron Flux for a 
Central Disturbance.. 
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